clear variables 

datapath = 'Output\';
filename = 'Chains.xlsx';
name = 'replication';

omega_bar = 0;

load(['parameters_22\','param_ini_',name,'.mat']);
load(['parameters_22\','param_end_',name,'.mat']);

region='R';

aux.n_m = 10^7;
aux.n_j = 50;

%% Set c=0 and resolve model
param_ini.c=0;
param=param_ini;
G_l=1+(param.c+param.K).*(param.r+param.zeta)./(param.omega_0);
G_u=(1-1./param.gamma_R_1)./(1-1./param.gamma_R_2).*(1+(param.c+param.K).*(param.r+param.zeta)./(param.omega_0));%param.s.*param.lambda+    
tmp_G= fzero(@ (G) (param.omega_0 + rho_fun(0,region,param).*(param.c+param.K))./param.omega_0.*phi_fun( G,region, param) ...
            -G.*phi_fun( G.^(-1), region, param), [G_l,G_u], optimset('display','off'));
param.G=tmp_G;

param.m_l = param.omega_0./(1-param.omega_1)./rho_fun(0,'R',param)./phi_fun(param.G,region, param);
param.m_e = (param.omega_0+rho_fun(0,'R',param).*(param.c+param.K))./(1-param.omega_1)./rho_fun(0,'R',param)./phi_fun(param.G.^(-1),region, param);
param.m_u = param.m_e;            
param.m_h = param.m_l;            
param.m_l_ini=param.m_l;
param.m_u_ini=param.m_u;
param.shock = 0.0001;

param_ini=param;


if omega_bar == 0
    % Omega_0 is acyclical
    param_end=param;
else
    % Otherwise to allow omega_0 to change use
    param_end.c=0;
    param=param_end;
    G_l=1+(param.c+param.K).*(param.r+param.zeta)./(param.omega_0);
    G_u=(1-1./param.gamma_R_1)./(1-1./param.gamma_R_2).*(1+(param.c+param.K).*(param.r+param.zeta)./(param.omega_0));%param.s.*param.lambda+    
    tmp_G= fzero(@ (G) (param.omega_0 + rho_fun(0,region,param).*(param.c+param.K))./param.omega_0.*phi_fun(G,param) ...
                -G.*phi_fun(G.^(-1),param), [G_l,G_u], optimset('display','off'));
    param.G=tmp_G;
    
    param.m_l = param.omega_0./(1-param.omega_1)./rho_fun(0,'R',param)./phi_fun(param.G,param);
    param.m_e = (param.omega_0+rho_fun(0,'R',param).*(param.c+param.K))./(1-param.omega_1)./rho_fun(0,'R',param)./phi_fun(param.G.^(-1),param);
    param.m_u = param.m_e;                      
    param.m_h = param.m_l;
end


R = (param.m_u/param.m_l)^(1/(2*(1-param.alpha)));
j=1:aux.n_j;
lambda = 1/8.*param.sigma^2.*(1./(1-param.alpha))^2 *(1 + (pi.*j./log(R)).^2);

tt = (1:100)';

dlnN_t = nan(size(tt));
for k = 1:length(tt)
    dlnN_t(k) = -1./(1-param.alpha).*(1-sum((2.*pi*j./((pi.*j).^2+(log(R)).^2)).^2.*(1-(-1).^j.*((R+R^(-1))./2)).*exp(-lambda.*tt(k))))*param_ini.shock;
end

Months = tt;
T = table([0;Months],[0;dlnN_t]);
writetable(T,[datapath,filename],'Sheet','Figure G','Range','A1')


%% Check

lnm = linspace(log(param.m_l),log(param.m_u),aux.n_m)';
phi_j   = exp( - 1/2*1/(1-param.alpha)*(lnm - log(param.m_l))).*(cos(pi*j.*(lnm-log(param.m_l))/(log(param.m_u)-log(param.m_l))) + (log(param.m_u)-log(param.m_l))./(2*pi*j*(1-param.alpha)).*sin( pi*j.*(lnm-log(param.m_l))/(log(param.m_u) - log(param.m_l))));
phi_t_j = exp( + 1/2*1/(1-param.alpha)*(lnm - log(param.m_l))).*(cos(pi*j.*(lnm-log(param.m_l))/(log(param.m_u)-log(param.m_l))) + (log(param.m_u)-log(param.m_l))./(2*pi*j*(1-param.alpha)).*sin( pi*j.*(lnm-log(param.m_l))/(log(param.m_u) - log(param.m_l))));
M_j = (phi_t_j'*phi_j*log(param.m_u/param.m_l))'./length(lnm);
coef_j = diag(M_j);
max(max(abs(M_j - diag(M_j).*eye(size(M_j)))))

int_adj = ((log(param.m_u)-log(param.m_l))/2*(1 + ((log(param.m_u)-log(param.m_l))./(2*pi*j*(1-param.alpha))).^2));
max(1-coef_j'./int_adj)

param.m_l_ini = param_ini.m_l;
param.m_u_ini = param_ini.m_u;
Q_ini=linspace(0,1,aux.n_m)';
Q_new = Q_ini_fun(Q_ini,param,aux);
f = [Q_new(2:end)-Q_new(1:end-1);Q_new(end)-Q_new(end-1)];
ff = f./sum(f)./log(param.m_u/param.m_l)*aux.n_m;


int_num = (phi_t_j'*ff*(log(param.m_u/param.m_l)))'./length(lnm);
int_analytic = (param.shock./(log(param.m_u/param.m_l))*(1-(-1).^j*exp( 1/2*1/(1 - param.alpha)*(log(param.m_u/param.m_l)))));
max(abs(int_analytic-int_num))

ff_new = ones(size(ff))./(log(param.m_u/param.m_l));
ff_new(1)=sum(ff_new((lnm>log(param.m_u)-param.shock)));
ff_new((lnm>log(param.m_u)-param.shock)) = 0;

int_num = (phi_t_j'*ff_new*(log(param.m_u/param.m_l)))'./length(lnm);
max(abs(1-int_num./(param.shock./(log(param.m_u/param.m_l))*(1-(-1).^j*exp( 1/2*1/(1 - param.alpha)*(log(param.m_u/param.m_l)))))))

